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Abstract. 

In the extreme violence of merger and mass accretion, compact objects like black 
holes and neutron stars are thought to launch some of the most luminous outbursts 
of electromagnetic and gravitational wave energy in the Universe. Modeling these 
systems realistically is a central problem in theoretical astrophysics, but has proven 
extremely challenging, requiring the development of numerical relativity codes that 
solve Einstein’s equations for the spacetime, coupled to the equations of general 
relativistic (ideal) magnetohydrodynamics (GRMHD) for the magnetized fluids. Over 
the past decade, the Illinois Numerical Relativity (ILNR) Group’s dynamical spacetime 
GRMHD code has proven itself as a robust and reliable tool for theoretical modeling 
of such GRMHD phenomena. However, the code was written “by experts and for 
experts” of the code, with a steep learning curve that would severely hinder community 
adoption if it were open-sourced. Here we present IllinoisGRMHD, which is an open- 
source, highly-extensible rewrite of the original closed-source GRMHD code of the 
ILNR Group. Reducing the learning curve was the primary focus of this rewrite, with 
the goal of facilitating community involvement in the code’s use and development, as 
well as the minimization of human effort in generating new science. IllinoisGRMHD 
also saves computer time, generating roundoff-precision identical output to the original 
code on adaptive-mesh grids, but nearly twice as fast at scales of hundreds to thousands 
of cores. 


PACS numbers: 04.25.D-, 04.30.Tv, 04.40.Dg, 07.05.Tp, 47.75.+f, 52.30.Cv, 95.75.Pq 
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1. Introduction &; Motivation 

Incident gravitational wave (GW) observations have the potential to address some of the 
most important unsolved problems in astronomy and theoretical astrophysics. These 
include testing GR in the strong field regime, determining the engine behind short-hard 
gamma-ray-bursts, constraining the equation of state above nuclear densities, revealing 
how compact binaries form and evolve, as well as uncovering the distributions of black 
hole spins and masses, just to name a few. However, it is well known that unless 
some coincident electromagnetic (EM) counterpart from the GW source is observed, 

GW interferometers alone may be unable to pinpoint the source position on the sky, 
hindering parameter estimation [5T1 H6l I6H [35] . Moreover, EM signals carry additional 
and complementary information about the source, lending potentially critical insights 
about the GW source and its environment. 

Thus detections of EM counterparts to GWs could be critically important in this 
age of “multimessenger” astronomy, and not solely when GWs are detected first. For 
example, it may be possible that an EM signal itself would imply a GW source, leading 
to targeted searches across the GW spectrum, from the nHz band in the case of dual 
AGNs, to the kHz band in the case of stellar-mass binaries and supernovae. Beyond 
coincident GW detections, EM transients linked to strong-field, dynamical spacetime 
phenomena may themselves greatly advance our understanding of black hole (BH) 
accretion phenomena and matter at extreme densities. 

However, without detailed theoretical models of EM counterparts to GW 
observations, our interpretation of observed EM counterparts may be severely limited. 
Constructing such models remains a central problem in theoretical astrophysics, for two 
key reasons. First, observable signals are often sensitive to fluid flows and gravitational 
fields spanning many orders of magnitude in lengthscale and timescale. Second, the 
equations governing the dynamics are highly complex and nonlinear, requiring the 
evolution of the full set of Einstein’s equations of general relativity, coupled to the 
equations of general relativistic magnetohydrodynamics (GRMHD). 

Thus numerical relativity codes capable of modeling multi-scale GRMHD flows 
promise to not only provide key insights into these important phenomena, but represent 
the starting point for more sophisticated modeling that includes advanced EM and 
neutrino radiation transport. More than a decade ago, the Illinois Numerical Relativity 
Group led by S. L. Shapiro was among the first to develop a dynamical spacetime 
GRMHD code [25] for uniform-resolution grids. Since then, this GRMHD code 
(henceforth, OrigGRMHD) has been significantly extended and improved. Its current 
version models multi-scale GRMHD flows via an adaptive-mesh-refinement (AMR) 
vector-potential formulation. By evolving the vector potential forward in time instead 
of the magnetic fields directly, this formulation guarantees the no-monopole constraint 
V • B = 0 is satisfied over the entire numerical grid, even when multi-scale magnetized 
fluid flows cross AMR grid boundaries. Notably, this formulation reduces to the 
standard, staggered Flux-constrained-transport (FluxCT) [TT] scheme on uniform- 
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resolution numerical grids [3D], "32] 13D] . 

OrigGRMHD’s reputation for generating models that address key unsolved 
problems in theoretical astrophysics has been built upon years of hard-fought 
development, as there exists no standard, proven algorithms for dynamical spacetime, 
multi-scale GRMHD modeling. Over the past decade, the code has been used to model 
a number of astrophysical scenarios, gleaning key new insights into these systems. For 
example, OrigGRMHD has produced state-of-the-art magnetized binary neutron stars 
(NSs) [4?, ,h7] and binary BH-NS |'2T. 3T[ [28', 33. [60] models. It was also used to simulate 
magnetized disk accretion onto binary BHs [371 m EH E], binary white dwarf-NS 
mergers [Ml 59] , magnetized, rotating NSs [29], magnetized Bondi accretion and 
magnetized hypermassive NSs [65] EZl 23, [68] [66, 38], just to name a few. The code was 
also recently extended, as a separate module, to solve the equations of GR force-free 
electrodynamics |EH and applied to model both binary black hole-neutron star [58] and 
pulsar magnetospheres in GR [62]. At each stage of its development, OrigGRMHD was 
subjected to a large battery of stringent test-bed problems 03 . which it had to pass 
before being used for applications. 

The field has matured considerably in the years since the first dynamical spacetime 
GRMHD codes were announced, and multiple groups now possess their own independent 
codes [25] [2] [39] [16] [30] [32] [63] [331 [50] [21] [55], most of which solve these equations on 
adaptive mesh refinement (AMR) grids. Given the time and effort required to extend 
such codes to model more physics, while still maintaining and improving the GRMHD 
modules, it seems clear that the community might benefit if more of us consolidated our 
efforts and adopted the same dynamic-spacetime GRMIID code. 

With its proven robustness and reliability in modeling some of the most extreme 
phenomena in the Universe, OrigGRMHD appears to be a good candidate for such 
community adoption if it were open-sourced. But despite its strong track record, 
OrigGRMHD was not written with community adoption in mind, instead being a code 
written “by experts and for experts” of the code. As such, the code lacked a number of 
features common to top-notch, widely-adopted open-source projects in computational 
astrophysics, including sufficient documentation and code comments, fine-grained 
modularity, a consistent coding style, and regular, enforced code maintenance (e.g., 
removal of unused and unmaintained features). 

Thus the OrigGRMHD core development team came to the understanding that 
unless these idiosyncrasies were fixed, open-sourcing the code would be unlikely to 
engender widespread community adoption. Thus in early 2013, it was decided to rewrite 
OrigGRMHD from the ground up, with a focus on the four core design principles of user- 
friendliness, robustness, modularity/extensibility, and performance/scalability. Slightly 
more than a year later, all of OrigGRMHD’s core algorithms had been rewritten and 
the new code, IllinoisGRMHD, was released. 

Just after the decision was made to rewrite OrigGRMHD in 2013, the first 
open-source, dynamical spacetime GRMHD code, called GRHydro, was released [50]. 
Originally forked from the dynamical spacetime, general relativistic hydrodynamics 
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(GRHD) Whisky code [Hj (not to be confused with its closed-source successor, 
WhiskyMHD [39]), GRHydro shares many of the same features of OrigGRMHD, 
including a number of reconstruction techniques. 

However, unlike OrigGRMHD/lllinoisGRMHD, GRHydro’s GRMHD scheme has 
not been developed to forbid the generation of monopoles (i.e., violations of the 
V ■ B = 0 constraint) when magnetized fluids flow across AMR grid boundaries. 

As accurate modeling of such multi-scale fluid flows is critically important in many 
astrophysical scenarios of interest to the community, GRHydro’s adoption by members 
of the community has been limited, primarily to those who simulate core collapse. 

Further, one of OrigGRMHD/lllinoisGRMHD’s key advantages is that these codes 
are capable of stably modeling GRMHD flows into black hole (BH) horizons over very 
long timescales, without the need for special algorithms that excise GRMHD data within 
the BH. By contrast, it seems that BH excision is an essential ingredient for stable 
GRMHD evolutions with GRHydro in the presence of black holes. To date GRHydro 
has been mostly used for core collapse (to a neutron star) simulations, in which no black 
hole is present. 

OrigGRMHD/lllinoisGRMHD have been demonstrated robust across a much wider 
range of long-term BH simulations, and manage to do so without excision. We conclude 
that making GRHydro’s GRMHD schemes as robust may require careful specification of 
boundary conditions on the excision surface coupled to an interpolation scheme across 
AMR level boundaries that respects the no-monopoles constraint, e.g. mm- 

IllinoisGRMHD was originally designed in a standalone sandbox to maximize 
portability to other parallel infrastructures, but currently adopts the latest Einstein 
Toolkit (ET)/Carpet AMR infrastructure. IllinoisGRMHD has been proposed for 
inclusion within the next ET release, and code review is underway. In the meantime, the 
ET community have graciously agreed to host the IllinoisGRMHD code, in anticipation 
of official incorporation upon completion of the code review processjj] 

In this paper, we present results from a number of code validation tests 
demonstrating that IllinoisGRMHD (1) produces results identical to OrigGRMHD, (2) 
possesses identical or significantly better scalability and performance than OrigGRMHD 
and GRHydro, (3) generates results in quantitative agreement with those of the 
GRHydro code, in the context of dynamical spacetime evolutions of Tolrnan- 
Oppenheimer-Volkoff (TOV) stars. 

The remainder of the paper is organized as follows. Section [2] outlines the 
formulation of the GRMHD equations solved by IllinoisGRMHD, Sec. [3] presents a 
basic overview of algorithms used within IllinoisGRMHD, Sec. [4] shows results from 
code validation tests, Sec. [5] demonstrates the outstanding performance and scalability 
of IllinoisGRMHD via benchmarks, and Sec. [6] summarizes results and describes future 
work. 

f Instructions for downloading, compiling, and using IllinoisGRMHD may be found here: 
http://math .wvu. edu/~zetienne/ILGRMHD/ 
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2. Basic Equations 


All equations presented below are in geometrized units where G = c = 1. In these units, 
Einstein’s equations become 

= 8v tT^, (1) 

where G^ v is the Einstein tensor and T^ v the total stress-energy tensor. IllinoisGRMHD 
solves the coupled Einstein-Maxwell equations assuming a perfect fluid stress-energy 
tensor for the matter and infinite conductivity (ideal MHD), by evolving via liigh- 
resolution-schock-capturing techniques the GRMHD quantities that comprise the stress- 
energy tensor T M ", acting as the source for Einstein’s equations. With these assumptions, 
the GRMHD evolution and constraint equations are derived from the following basic 
equations: 


(i) Conservation of baryon number 

V^po^) = 0, 


( 2 ) 


where is the covariant derivative associated with the spacetime metric tensor 
Po is the fluid rest-mass density and is the fluid four-velocity. 

(ii) Conservation of energy-momentum 




(IT = T 

U i -L III/ -L , 




+ T, 


flU 


(3) 


where is the sum of the perfect fluid T^ atter and electromagnetic stress-energy 
tensors T™ in the ideal MHD limit (u^F^ v = 0). 

(iii) Homogeneous Maxwell’s equations 

= ~^= d v (^F*n = 0, (4) 

where F^ the Faraday tensor, F* ,lL ' = (1/2 )e tJ ' upa F prr its dual (e Mi/pCT is the Levi- 
Civita tensor), and g the determinant of g /w . 


As written, Eqs. (| 2 |), (13|), and ()4]) for the plasma, as well as Eq. ([T]) for the spacetime 
metric, are not particularly well-suited for numerical evolutions, so we choose special 
formulations of them. For the spacetime metric evolution, we choose an initial value 
formulation built upon first splitting the 4-metric g^ u into the standard 3+1 Arnowitt- 
Deser-Misner (ADM) form [5]: 

ds 2 = g^dx^dx 1 * = — a 2 dt 2 + 7 ij(dx l + f3 l dt)(dx J + ftdt). (5) 


Here, a is the lapse function, /+ the shift vector, and 7 ^ the three-metric on spacelike 
hypersurfaces of constant time t. This basic decomposition of the 4-metric can be used 
to split the Einstein equations (JT]) into a set of evolution equations and a set of constraint 
equations that the dynamical variables must satisfy for all times—similar to Maxwell’s 
equations—with projections of along and normal to the 3D spatial hypersurface 
existing as source terms (see, e.g., na for a detailed discussion and references). This 
original formulation of the Einstein equations is known as the ADM 3+1 formulation 



IllinoisGRMHD: An Open-Source, User-Friendly GRMHD Code for Dynamical Spacetimes 6 


of GR. A number of 3+1 formulations can be derived from the ADM formulation 
and are useful for solving the Cauchy problem for the Einstein equations. For the 
purposes of this paper, we choose the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) 
formulation mm, which introduces an auxiliary dynamical variable and conformal 
scalings for the dynamical variables, casting the evolution equations in a form that 
allows for stable, long-term and accurate numerical integration. 

To update from one time slice to the next in a simulation, IllinoisGRMHD solves 
Eqs. 02 ]), ©, and the spatial component of Eq. (J4]) in the ideal MHD limit (u l _ l F fll/ = 0), 
as written in conservative form (see e.g. [25]): 


d t C + V • F = S , ( 6 ) 

where F is the flux vector, C = {/?*, f, Si, B 1 } the vector of conservative variables, and 
S the vector of source terms. These vectors depend directly on the “primitive” variables 
P = {p 0 , P,v\ B 1 }, where P is the pressure, v l = u l /u° the fluid three-velocity, and B l 
are the spatial components of the magnetic field (R M ) measured by normal (or Eulerian) 
observers with four velocity = (1, —/3 l )/a (and is normal to the spatial hypersurface, 
B^n^ = 0, as well as the metric and its derivatives. In particular, C may be written in 
terms of P , a and the metric as follows: 


' p * ’ 


Oiy/lPoU 0 

f 


cEyyT 00 - p* 

1 

*cq 

_ 1 


(p*h + au° y/^fb 2 )ui — a.y/^b°bi 
VlB\ 


( 7 ) 


where 7 is the determinant of the 3-metric 7 ^, h — 1 + e + P/ p 0 is the specific enthalpy, 
with e the specific internal energy, and U = B^/\fpf with the magnetic field field 
measured by an observer comoving with the fluid. Here the total stress-energy tensor 
T can be written as follows (see [HI [25] [24] for further details and derivations) 


T M " = (p 0 h + b 2 )vPu v + 



g^ - UP. 


( 8 ) 


Our choice of fluid 3-velocity v l = u l /u° as a primitive variable is consistent with 
a number of GRMHD codes, such as [5311261111]. However, our v l differs from that of 
other GRMHD codes (e.g., [501 E, E, 39], just to name a few) that have adopted the 
Valencia formalism C 2 li, which adopts the fluid 3-velocity v ^ as measured by normal 
observers (also referred to as the Eulerian 3-velocity), defined as: 


u a = au°(n a + u“ n )) 


u 


V H = 


,0 


+ 




( 9 ) 

( 10 ) 


au u a 

Note that v 1 ^ is orthogonal to the normal vector to the 3D spatial hypersurface 


v 


= 0. In terms of the fluid 3-velocity used by IllinoisGRMHD v l , v^ can be 


(n) 

written as 


+ = - 7 +n 

v a 


( 11 ) 







IllinoisGRMHD: An Open-Source, User-Friendly GRMHD Code for Dynamical Spacetimes 7 

Note that this difference in the 3-velocity variable may account for some of the 
differences in numerical results observed between IllinoisGRMHD and (the Valencia- 
based) GRHydro in Sec. QJ as Valencia-based codes reconstruct v 1 ^ instead of v l . 

Writing the GRMHD evolution equations in conservative form offers a number of 
numerical advantages. First, without the source terms S, it guarantees conservation 
of total rest mass (f v p*d 3 x), energy ( f v fd 3 x ), and momentum (f v Sid 3 x) to roundoff 
error. When the source terms are accounted for, total ADM mass and momentum 
are conserved to within truncation error. Second, it enables us to attach easily an 
approximate Riemann solver, yielding a state-of-the-art high-resolution shock-capturing 
(HRSC) numerical scheme, designed in part to minimize oscillations near shocks. Such 
oscillations are generated by approximating the flux derivative across shocks by smooth 
functions. Third, the conservative form, coupled to the HRSC scheme guarantees the 
shock jump conditions (Rankine-Hugoniot) are satisfied. From an empirical perspective, 
finite-volume conservative formulations, as implemented in IllinoisGRMHD, have been 
shown superior at handling ultrarelativistic flows when compared to advanced artificial 
viscosity (AV) schemes [JU 3], despite the fact that AV schemes are typically superior 
in terms of ease of implementation and computational efficiency. 

A key ingredient in a robust GRMHD code is the proper treatment of the magnetic 
induction equation, which is derived from the spatial components of Eq. (J4|). In 
conservative form, these may be written: 

d t B i + dj (ViT - v l B j ) = 0. (12) 

If the induction equation is directly evaluated in a numerical code without special 
techniques, numerical truncation errors that violate the divergence-free or “no¬ 
monopoles” constraint 

ditf = 0 (13) 

will be generated. Note that this equation is simply the time component of Eq. (j4j). 
Maintaining satisfaction of this constraint as the magnetic fields are evolved forward 
in time [through direct evaluation of Eq. (fT2l) ] happens to be a nontrivial endeavor, 
particularly on AMR grids. Our solution |3_Q, E2j is to evolve the magnetic 4-vector 
potential A p instead of the magnetic fields directly, so that 

A p = + Ap, and (14) 

B l = U^djAk, (15) 

where A M is purely spatial (A^n^ = 0) and $ is the EM scalar potential. Here, e l]k is 
the standard permutation symbol, equal to 1 (-1) if ijk are an even (odd) permutation 
of 123, and 0 if one or more indices are identical. Special finite difference operators for 
the vector potential are defined in IllinoisGRMHD so that the divergence of a curl is 
zero to roundoff error, which implies that the divergence of Eq. (TT5h is zero and Eq. (fT3l) 
is satisfied automatically, even on AMR grids. 

In terms of Ai, the induction equation (fT2l) becomes 

d t Ai - e ijk v 0 B k - 'Vo'!’ - 0 3 A,). 


( 16 ) 
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What remains is to choose an EM gauge, and IllinoisGRMHD chooses the “generalized 
Lorenz gauge condition” by default that was introduced by the Illinois Relativity group 
in [36]. The covariant version of the condition is where £ is a parameter 

with dimensions 1/Length, chosen carefully so that the Courant-Friedrichs-Lewy (CFL) 
factor remains satisfied. Typically £ is set to 1.5/Ai max , where At max is the timestep of 
the coarsest refinement level. This gauge choice therefore yields the additional evolution 
equation 

dt[y/l$] + djia^SyA 3 - /^[VtL]) = -£a^/ 7 $. (17) 

Note that IllinoisGRMHD evolves not $ but yWL as the EM gauge variable. 

With the exception of this purely gauge evolution equation and the vector-potential 
induction equation, all other GRMHD evolution equations are written in conservative 
form and are solved via a HRSC scheme, as described in Sec. [31 For completeness, 
the remaining set of evolution equations evolved by IllinoisGRMHD are written in 
conservative form [Eq. (JHJ)] as follows 



P * 


p*v j 


0 

d t 

1 

ICO 

+ 

a 2 v / 7T 0j ' - p+v 3 
a^T 3 i 

— 

s 

_ _ 


where 


s = a^7 (' T “W* + 2T Ul /3 J + T l3 )K i:j - (T uu /T + T^a 


-rOi 


00 I 


-1 On 


(19) 


and K t j = — £ n ^ij/2 is the extrinsic curvature, where £ n designates the Lie derivative 
along the hypersurface normal vector n (see e.g. ra for more details). 

Finally, to close the system of equations, the EOS of the matter must be specihed. 
IllinoisGRMHD currently implements a hybrid EOS of the form [ 43] 


P(p0; ^) Pcold(Pd) T (Tth l)p 0 [e Oold(Po)] > 


( 20 ) 


where P co id and e co id denote the cold component of P and e respectively, and T th is a 
constant parameter which determines the conversion efficiency of kinetic to thermal 
energy at shocks. The function e co id(po) is related to P co id(Po) by the first law of 
thermodynamics, 

tc„,d(Po) = [ PPfP-dp, . (21) 

P 0 

All functions within IllinoisGRMHD support piecewise-defined P co id(po) (the so-called 
“piecewise polytrope” EOS) with up to nine different polytropic indices, except for the 
conservatives-to-primitives solver, which currently supports only one. 

In all code tests presented in this paper, the T-law EOS P = (T — l)po e is adopted. 
This corresponds to setting P co i d = (T — l)poe co id in Eq. (120]) . which is equivalent to 
Pcoid = O°o ( w ith constant /?), and T t h = T. In the absence of shocks and in the initial 
data used for our tests, e = e co id and P = P co id- 
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2.1. Outer Boundary Conditions 

We apply outer boundary conditions to primitive variables po, P , and v\ which enforce 
a zero-derivative, “copy” boundary condition of these quantities at the outer boundary, 
except when this results in a positive incoming velocity from the outer boundary. Hence 
we refer to these as “outflow” boundary conditions. As our outer boundary exists as 
a rectangular box, if for example v x < 0 at a given point on the x = x max boundary 
plane after applying the zero-derivative boundary condition, we set v x = 0 at that point. 
Similarly, at a given point on the x = x min boundary plane, if v x > 0 after applying the 
flat outer boundary condition, v x is set to zero. The same strategy is applied to the 
velocities for all other outer boundary faces. 

We also apply outer boundary conditions to A^, linearly extrapolating values 
to the outer boundary. To avoid problems caused by reflections of A M waves from 
these imperfect outer boundary conditions, the unit-bearing (i.e., not dimensionless) £ 
parameter in the Ai evolution equation is set to some nonzero, positive value, typically 
1.5/Ai max , where At max is the timestep of the coarsest refinement level. 

3. Basic Algorithms 

Since its initial stages, one of the primary objectives driving the development of 
IllinoisGRMHD has been to remove nonrobust algorithms and obsolete code from 
the original GRMHD code of the Illinois group, resulting in a reliable state-of-the- 
art piece of software that is more compact and easier for beginners to learn and 
extend. To this end, all of the obsolete code and functionality proven to be nonrobust 
in typical dynamical spacetime evolutions has been stripped from the code, keeping 
within IllinoisGRMHD only the set of algorithms used in all of the Illinois group’s latest 
GRMHD publications. Further, the core algorithms themselves have been rewritten 
into a uniform coding standard, with large amounts of duplicated functionality replaced 
with a small, optimized library of functions. This section reviews the basic algorithms 
that comprise IllinoisGRMHD. 

Here the algorithms that comprise the basic components of IllinoisGRMHD are 
introduced, in the order in which they are called. At the beginning of the first timestep, 
the variables {P, po,v l , B l , A^, $} must be defined at every gridpoint. The following 
outlines the basic steps in which these variables are updated at all gridpoints, in 
preparation for the next timestep. All updates are performed by IllinoisGRMHD unless 
otherwise specified. 

(i) First, the flux and RHS terms in Eq. (fTTl) and Eqs. (fT8j) . for the set of evolution 
variables E = {p*, Si, t, A it [^/ 7 $]}, are evaluated. Three separate algorithms are 
employed in this step: 

(a) A HRSC evolution scheme is used to compute the flux terms of the p*, S tl and 
f evolution equations, as defined in Eqs. fjl8jl ). This scheme, as well as the 
technique used to compute the source terms related to spacetime curvature in 
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these equations, is described in Sec. 13.11 

(b) Unlike the other primitive variables, Ai and B l are defined on staggered 
gridpoints. Further, our evolution scheme is constructed to produce identical 
output to the standard, staggered constrained transport scheme of [34] . As 
detailed in Sec. 13.21 this makes the HRSC scheme for updating Ai a bit 
more involved than the HRSC scheme for evolving the unstaggered densitized 
density, momentum, and energy variables. 

(c) The evolution of the (staggered) EM gauge quantity [y / 7 < h] is not based on 
a HRSC scheme, as this quantity generally does not exhibit sharp features in 
our simulations. Its evolution algorithm is summarized in Sec. 13.31 

(ii) Time derivative data from all evolved GRMIID variables are then passed to the MoL 
(Method of Lines) thorn, which iteratively integrates the evolved variables forward 
in time. MoL is capable of managing a number of iterative, explicit time integration 
techniques, of which we typically choose the four-iteration Runge-Kutta fourth- 
order (RK4) scheme, both in IllinoisGRMHD and the chosen spacetime evolution 
thorn. 

(iii) Evaluating the time derivatives of all evolved GRMHD variables requires three 
ghostzones at the outer boundary of each AMR grid. The ghostzones at the 
outermost boundary are filled at each RK4 iteration, using the outer boundary 
update procedure outlined in the next steps. However, ghostzones at each internal 
AMR grid boundary are allowed to accumulate until the end of the fourth RK4 
iteration. Since RK4 consists of four iterations, this yields a total of 3 x 4 = 12 
AMR grid boundary ghostzones that must be filled at the end of each full RK4 
timestep. To fill these 12 ghostzones for all evolved variables at the end of the 
fourth RK4 iteration, prolongation and restriction operators are applied, which 
interpolate between different levels of refinement in both space and time. Third- 
order, line-averaged Lagrange prolongation/restriction is performed on Ai and a/mL, 
and fifth-order Lagrange prolongation/restriction is performed on all other GRMHD 
evolved variables. 

(iv) Next, linear-extrapolation outer boundary conditions are applied to Ai and 

as described in Sec. 13.51 Then B l is computed from Ai at all gridpoints, as described 
in Sec. 13.21 

(v) The conservative variables p*, Si, f, and B l have at this point been updated at all 
needed gridpoints, except at the outer boundary. However, the primitive variables 
{P, po, v\ B *} do not yet exist at any gridpoints. As these variables are required at 
the outset of the next iteration, a conservative-to-primitives solver is called next, 
which at its heart employs a Newton-Raphson-based root-finder to invert Eqs. (171). 
computing primitive variables at each gridpoint based on the conservative variables 
at that point. Additionally, there are a number of consistency checks applied both 
before and after this solver is called. The procedure is outlined in Sec. 13.41 

(vi) Next, zero-slope, outflow outer boundary conditions are applied to the set of 
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primitive variables {P, p 0 ,v 1 }, as described in Sec. 13.51 After this step, all variables 
{P, po, v\ B\ A needed to repeat this process have been defined at all gridpoints, 
so to proceed to the next timestep or RK4 iteration, we simply loop to (i)(a). 

Values for the conservatives at the outer boundary are not strictly required for the 
evolution, but are set anyway, based on the primitive variables, in case a diagnostics 
utility might require that conservatives be set at the outer boundary. 

To conclude this introduction to IllinoisGRMHD’s basic algorithms, Sec. 13.61 
describes how IllinoisGRMHD connects to the rest of the ET and its spacetime metric 
evolution modules. 

3.1. Evolution of p* ; S t , and t 

To evolve the GRMHD variables, IllinoisGRMHD first evaluates the source terms of the 
time evolution equations for {/?*, S',, r} in Eqs. (fl8|) . Derivatives of the spacetime fields 
appear in the source terms, which are evaluated via standard, second-order (default) or 
fourth-order finite differences. Next, the fluxes are computed via a second-order finite- 
volume HRSC scheme. Since point values of gridfunctions and their volume averages 
are the same to second order, our hnite-volume scheme will converge at the same order 
as a finite-difference scheme. We choose a hnite-volume scheme, as it will enable us 
to more rapidly move to a higher-order method in future releases of IllinoisGRMHD, 
following the strategy outlined in [30]. 

Computation of the flux term V • F = d m F m in a given direction i G { x , y, z} is 
performed with our second-order hnite-volume scheme in two steps, as detailed below. 

First, the Reconstruction Step computes values for the primitive variables at cell 
interfaces (between gridpoints) along direction i. Then the Riemann Solver solves the 
Riemann problem via an inexpensive, approximate algorithm, ensuring the conservative 
variable huxes between gridpoints are appropriately constructed along direction i, even 
in the presence of discontinuities or shocks. Upon completing the Riemann solver step 
for a given hux direction i , the process is repeated in the other two directions until 
d m F m has been evaluated and summed in all three spatial dimensions {x,y,z}. 

The Reconstruction Step for p*, 5), and f: IllinoisGRMHD employs the 
Piecewise Parabolic Method (PPM) [T7] . incorporating the original battening and 
steepening procedures to reconstruct P at the right (Pr) and left (Pl) sides of each 
grid zone interface, along direction i G {x,y,z}. The version of PPM used within 
IllinoisGRMHD is designed to maintain third-order accuracy, except at discontinuities 
or shocks and at local minima and maxima. As in the GRHydro code [(50] the flattening 
procedure within PPM was simplified to decrease the number of required ghostzones 
within PPM from four to three. 

After PPM reconstruction evaluates Pr,l along a given direction i and the metric 
values have been interpolated to each grid zone interface at fourth-order (default) or 
second-order accuracy, the huxes Fr,l are then immediately evaluated via (J7]) and 
Eqs. (HU). 
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Next, for appropriate handling of fluxes across a given cell interface, the Riemann 
problem must be solved. 

The Riemann Solver for p*, Si, and f: 

The first step in solving the Riemann problem along direction i G {x, y , z} is to 
compute the maximum (+) and minimum (—) characteristic speeds c± L at each cell 
interface, approximating the general GRMHD dispersion relation (Eq. 27 of [38]) with 
the following, simpler expression: 


(jJ nr r, 


VA + Cs(l - V\) 


ki 


( 22 ) 


Here, cu cm = —/qpG is the frequency and the wavenumber of an MHD 

wave mode in the frame comoving with the fluid, where K tl is defined as the projection 
of the wave vector k v onto the direction normal to u v : K M = ( g^ + Ufj,u u )k u . c s is the 
sound speed, and v\ is the Alfven speed, given by 


a i = 
0-2 = 
03 = 
„2 / 


VA V Poh + b 2 ' (23) 

With these dehnitions, the approximate dispersion relation [Eq. Q22]) ] may then be 
solved along direction i, noting that the wave vector along this direction in the comoving 
frame is given by k M = (—co, kjSj) and the wave (phase) velocity by c± = uj/(kj6j). The 
dispersion relation can then be written as a quadratic equation for c±: aic± + a 2 C± + a 3 = 
0, with ai given by 

(l-vl)(uy-vlg™, 

2v$g i0 - 2u i u°{l - vl), (24) 

(! ~ vl)^) 2 - vlg™, 
and vl = v\ + Cg(l — vjf). 

Though it makes c± simple to compute, this dispersion relation overestimates the 
maximum characteristic speeds by a factor < 2, which has the net effect of making the 
code more dissipative. Though additional dissipation may smear important physical 
features in our GRMHD flows, it also acts to help stabilize evolutions. Note that this 
approximate dispersion relation is widely used in multiple codes within the GRMHD 
community (e.g., WhiskyMHD [[39], GRHydro [50], HARM3D [53]). 

Once the maximum and minimum speeds c± have been computed at left and right 
faces, the standard Harten-Lax-van Leer (HLL), approximate Riemann solver 
then applied to compute fluxes for the three conservative variables U = {p*, f, Si}: 

c~F r + c + F l - c + c~(U R - U L ) 


is 


^HLL _ 


(25) 


where cA = ± max(0, c^, c^), and U RtR are the conservative variables U computed from 
the right and left reconstructed primitive values Pr,l, respectively. 

Upon computing the HLL flux at cell interfaces, the final step in evaluating the 
flux terms in the evolution equations of p*, S t , and f [Eqs. (|T8j) ] is to differentiate the 
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computed HLL flux terms along the same direction in which they were evaluated. After 
computing the HLL flux in the x-direction we calculate the x-derivative of the flux as: 

(d x F x ) i>jik = ( 26 ) 

The remaining y and 2 -terms in the d m F m sum are added to the sum as reconstruction 
proceeds along the y and 2 -directions, respectively. 

As the source terms of Eqs. (jT8|) have already been computed, to complete the 
evaluation of dtp*, d t Si, and d t f, all components of the d m F m sum are then subtracted 
from the source terms. These data are then passed to the MoL (Method of Lines) thorn, 
which is capable of managing a number of explicit time-stepping techniques. Although 
MoL supports a total-variation diminishing third-order Runge-Kutta time integrator, we 
typically choose the Runge-Kutta fourth-order (RK4) scheme for all evolution variables, 
both in IllinoisGRMHD and the chosen spacetime evolution thorn, as we find that RK4 
minimizes the total error when evolving both the spacetime and the fluid. 

In parallel with evaluating the flux and source terms for p*, Si, and f, 
IllinoisGRMHD employs a vector-potential-based constrained transport scheme to 
evolve the magnetic fields, which is detailed in the next section. 

3.2. Vector-Potential-Based Constrained Transport Scheme 

Constrained-transport schemes maintain V • B = 0 through careful finite differencing of 
the magnetic induction equation flux terms [Eq. (1T21) ]. Such schemes have proven highly 
robust in the context of strongly-curved spacetimes; in particular those inhabited by at 
least one black hole. These schemes are most commonly and most directly applied 
in the context of uniform-resolution grids. However, their use with AMR grids can 
be complicated, as maintaining the divergenceless constraint at refinement boundaries 
requires that special interpolations be performed during prolongation/restriction. Such 
prolongation/restriction operators have been devised [H|T0j, but must be fine-tuned to 
the particular AMR implementation. 

IllinoisGRMHD applies an alternative constrained-transport scheme, introduced 
by [2D]. In this scheme, the magnetic induction equation (fl2l) is recast as an evolution 
equation for the magnetic vector potential [Eq. (fTfil) ]. This scheme has two important 
advantages. First, it produces identical output to the standard, staggered constrained- 
transport scheme on uniform resolution grids and thus shares its robustness. Second, 
evolving the vector potential enables us to use any interpolation scheme at AMR 
refinement boundaries without introducing nonzero divergence to the magnetic fields, 
so long as we compute B l from the interpolated A{. 

The remainder of this section details the staggered constrained-transport scheme 
adopted within IllinoisGRMHD. First, we define the staggerings of individual 
gridfunctions and the computation of B rn from A rn . Then the technique of 
reconstruction [Eq. (TT2l) ] on staggered cell faces is outlined, and finally the Riemann 
solver is described. 
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Computation of B m from A m : 

In employing the standard, staggered constrained transport scheme, magnetic fields 
are defined at gridpoints that are staggered with respect to other conservative variables, 
as specified in Table [0 Notice that A m is staggered so that B m may be computed 
immediately from Eq. (fT5l) using the following finite difference representation, accurate 
to second order: 


TJX 

i+\,j,k 


V7i + ij t k(®yAz) 


where 




i+Uj,k 


= exp 


i+Uj,k 




V+o,j,k 


6 x g A 


(d y A z ) 




AA l+lt - A 


and 


(.dzA 


\i+i,j,k _ 


Ay 

i+\,3,k+\ 

-Gi-v 


A 


Az 


(27) 


(28) 


(29) 


(30) 


Here, = (1/12) logy is the BSSN conformal exponent. Using Eq. (1271) as a template, 
B y and B z can be written via straightforward permutation of vector indices {x,y,z}, 
accounting for the appropriate staggerings. Our finite differencing scheme is specified 
so that the divergence of a curl is identically zero to roundoff error. In this way, 
B m is guaranteed to be divergenceless at all but the outermost ghost-zones on any 
given refinement level, so long as A m is computed at all points. As with the other 
conservative variables, reconstruction of the flux terms for A m requires three ghostzones 
(as discussed in the next section), so the prerequisite step of computing B m from A m 
adds an additional ghostzone, bringing the total number to four. However, we have found 
that application of a copy boundary condition on B m to the outermost gridpoint on each 
refinement level, coupled to the use of only three ghostzones, results in qualitatively 
identical results to runs that use four ghostzones. We find this to be the case even 
in the most stringent tests, such as a magnetized BH accretion disk crossing multiple 
refinement boundaries, as in mmm- Thus by default, we have used 3 ghostzones 
in all GRMHD simulations. 

Flux Reconstruction of the Induction Equation: 

Accounting for staggerings, the evolution equation for A z (dropping the EM gauge 
terms to focus on the flux term of the induction equation) is given by 


O t A z - -b i+ y + T k , 


(31) 


where 


£ z = — v x B y + v y B 2 


(32) 


is the flux term in the standard magnetic induction equation ([T2]h Following [20], 
we compute this flux term to staggered cell faces for A z and then evaluate the HLL 
flux generalized for staggered grids. As £ z does not appear within a derivative of 
the A z induction equation (as it does in the B z induction equation), the flux is not 
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Table 1 . Storage location on grid of the magnetic field B l and vector potential Ay. 
Note that P is the vector of primitive variables {po,P,v 1 } 


Variable(s) 

storage location 

Metric terms, P, p*, S), f 


B x , B x 

(i + i,j,k) 

B y , B y 

(i,j + b k ) 

B z , B z 

(i,j,k+ \) 

A x 

(b j J r\,k J r |) 

Ay 

(® + \i3i k + |) 

A z 

(* + \,3 + b k ) 


(i + |)3 + ^ k + \) 


directly finite-differenced prior to passing the right-hand side of d t A z to the time- 
stepping routines. Instead, the spatial finite difference is computed after each RK4 
iteration when B m is computed from A m . Critically, the order in which spatial and 
temporal derivatives are evaluated is the only difference between the standard, staggered 
constrained-transport scheme and our vector-potential based staggered constrained- 
transport scheme. And since spatial and temporal derivative operators commute within 
IllinoisGRMHD’s current framework, both schemes are identical on uniform meshes. 

Returning to the evaluation of £ z , recall that primitives such as v l are defined at 
grid points (■ i,j , k), so computing the value £ z at {i + \,j-\-\,k) requires two successive 
one-dimensional reconstructions of v x and v y : first in the x or ^-direction and then in 
the y or x-direction, respectively. B x and B y already exist on staggered gridpoints (see 
Table [Tj) , requiring only a single reconstruction in the y and x direction, respectively. 
Reconstruction is handled via the same PPM scheme as described in Sec. 13.11 

Approximate Riemann Solver for Ap. The standard HLL formula (1251) for £ z , 
generalized to the appropriate staggered gridfunctions, is given by: 




2\HLL 


C t C y £ LL + C x C y £ LR + C„ £ fiy, + C r C„ £f 


-'x 0 RL 


-x L y °RR 


(' C x + C x ) ( C t + C y ) 


+ 


C x + C x 


-(B y R -B y )-^L-(B x -Bf) 


C+ 4- c~ 
y ' u y 


(33) 


In the above formula, ££ R denotes the reconstructed left state in the x-direction and 
right state in the ^/-direction. Other symbols involving £ z are interpreted in the similar 
fashion. B R (B\f) denotes the reconstructed right (left) state of B y in the x-direction, 
and B r (Bf) denotes the reconstructed right (left) state in the y-direction. The cff 
and c± should be computed by taking the maximum characteristic speed among the 
four reconstructed states. However, we set them equal to the maximum over the two 
neighboring interface values for simplicity, as suggested in [20], using the technique 
described in Sec. 13.11 to estimate the speeds. The formula for (£ X ) HLL is obtained from 
Eq. ()33j) by permuting the indices z —> x, x — > y and y —>■ z, whereas the formula for 
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(£ y ) HLL is obtained from Eq. (Hoi) by permuting the indices z —>■ y, x —$■ z and y —> x. 


3.3. Evolution of the densitized EM scalar potential [■ x /7 < h] 


Incorporating the staggering of the EM gauge variable [■ v /7 < h] (as specified in Table [Q), 
the evolution equation for [y / 7 < h] (Eq. [T7]) may be written: 


dt[Vl*]i+y+hMh = 


-d m {a^ mn An) + d m (/S m [v/7$]) - (34) 

Term (1) 


Term (2) 


Term (3) 


where ip is the standard BSSN conformal factor and the relations r y mn = -0 _4 y mri 
and y/j = ip 6 have been applied. The left-hand side of the equation is evaluated at 
(i + \,j + |, k + |), yet A m , [y^ 7 $], and metric quantities on the right-hand side of this 
equation all possess different staggerings. Thus special care must be taken so that the 
derivatives on the right-hand side (RHS) of this equation are evaluated at gridpoints 
[i + j + k + |). To accomplish this, quantities within the RHS derivatives are hrst 
interpolated to consistent points prior to evaluation of the derivatives. Note that this is 
strategy differs from the evolution of other GRMHD variables, in that no reconstruction 
is applied. Methods for computing these terms on the RHS are as follows: 

Term (1): For the ^-derivative, all quantities within the derivative operator 
([aip 2 ]^f mx A m ) are hrst interpolated to (i,j + + |). At second-order accuracy, 

interpolations to staggered gridpoints are trivial, requiring only averages of neighboring 
unstaggered points. For example, interpolation of aip 2 from ( i,j,k ) to (■ i,j + + |) 

is performed in two steps: 

(i) , k = |(M’ 2 kyfc + [oap 2 \i, j+ 1 ,k) 

(ii) [uip\ j+ ±, k+ l = + [onp 2 ] ijj+ i k+ 1) 

Once [aip 2 ], 7 m R A y , and A~ have been interpolated in this way to (■ i,j + |, k + |), the 
derivative d x ([aip 2 ]*/ xm A m ) is computed to second order as follows 


d x ([aiP 2 ]r m A r , 


(N 2 ]f m 4) 


i+lj + Rfc+a 


([aip 2 ]^ xrn A r 




' , *+2>t+2’ fc +2 /\x 

Other derivatives in the sum d m (aip 2 'p rnn A n ) are computed in the same fashion. 

Term (2): The computation of this term is made easier by the fact that [^/7$] 
is staggered at {i + |, j + k + |) already. So to evaluate the derivative, /3 m is hrst 
interpolated from ( i,j,k ) to (i + |, j + k + |) using the same interpolation strategy 
as with Term (1). Next, notice that this term is basically a shift advection term on 
the EM gauge quantity Such advection terms are typically upwinded within the 

metric evolution thorn, so for consistency we apply the same upwinding strategy when 
evaluating this derivative: 


dml 8 m [v/ 7 *]) 


d + (/?’"[ v /t»]) 


if /3 m < 0 , 
otherwise 


(36) 
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where the second-order operators are 


(D x f) j_|_I 


(37) 

2Ax 

f)i+^ ,j + i,fc+i 

~fi+ fc+i + 4/ i+ | ii+ i ifc+ i - Sf i+ i J+ i tk+ i 

(38) 

2Ax 


for the derivative in the ^-direction. Derivatives in the y and ^-directions follow in a 
straightforward fashion. 

Term (3): The computation of this term is also made easier by the fact that [y'qkl?] 
is staggered at + + + already. So to evaluate it, only a must be interpolated 
from ( i,j , k) to (i + 4, j + k + |) using the same interpolation strategy as with Term 
( 1 )- 

3-4- Conservatives-to-Primitives Solver 

After the conservative GRMHD variables have been updated at all gridpoints, with 
boundary conditions and prolongation/restriction operators applied, the primitive 
variables must then be computed from the conservative variables. This is not a trivial 
endeavor, as the conservative variables generally depend on the primitive variables in 
a nonlinear way, requiring the implementation of a root-finding method. To this end, 
IllinoisGRMHD employs the two-dimensional Newton-Raphson solver of [52l 53], 

Truncation errors originating from spatial and temporal finite differencing, as well 
as interpolation, prolongation, and restriction operations can push the evolved GRMHD 
quantities to unphysical values, resulting in either unphysical values for the primitive 
variables or no values at all. For definitions of unphysical values of the GRMHD 
quantities please see the Appendix A of [28j. So prior to calling the two-dimensional 
Newton-Raphson solver, we perform a number of checks that determine whether the 
conservative variables are in a physically valid range. If they are not, they are modified 
prior to calling the root-finder. Even with these checks, the Newton-Raphson solver 
will occasionally fail to find a root. This is very rare, and almost always occurs in a 
low-density atmosphere or inside a black hole. In such an instance, we set the pressure 
to P co id, which guarantees a successful inversion. The implementation of these checks 
and modifications have been described in detail in Appendix A of [28]. 

After the Newton-Raphson solver has successfully found a set of primitives, the 
primitives are checked for physicality, and if they are not in the physical range, they 
are minimally modified until they return to the physical range. First, if the velocity 
is found to be superluminal, the speed is reduced to IllinoisGRMHD’s default Lorentz 
factor limit, which is set to lb = 10, where W is the Lorentz factor of the fluid as 
measured by a normal observer. Next, IllinoisGRMHD does not include any cooling 
mechanism, which means that for evolutions adopting a T-law equations of state, the 
pressure should not physically drop below P co id- So a pressure floor of 0.9P co id is imposed. 
Increasing this floor to P co id exactly results in large central density drifts in TOY star 
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evolutions. Simulations can crash in the other extreme, if P/P co id becomes too large. 
This typically only happens in very low density regions or inside black holes. So at 
densities po < 100p atm or deep inside black hole horizons, a ceiling on P of 100P co id is 
enforced (see Appendix A of [28] for more details). 


3.5. Outer Boundary Conditions for A i} [ 1 / 7 $], P, Po, and v 1 


Updating evolved variables within IllinoisGRMHD requires three ghostzones per RK4 
iteration, and at the end of each iteration, outer boundary conditions are applied to A;, 
[.^/ 7 $], P , po, and v l fill these ghostzones. The algorithm applies the outer boundary 
conditions in all directions, from the innermost gridpoint outward, as follows. For 
example, in the positive x— direction, the first outer gridpoint * + 1 is defined as 


f Ei, if E e {P, po, v y , v z }, or E = v x and v x > 0 

E i+ 1 = < 0, if E = v x , and v x < 0 (39) 

( 2 Ei — Ei_i, if E G {["^/7<h], A x , A y , A z } 


And for the negative x— direction, the first outer gridpoint z — 1 is defined as 


Ei~ 1 — < 


Ei, 

0, 

2 E; 


E 


i+l) 


if E 6 {P, po, v y , v z }, or E = v x and v x < 0 
if E = v x , and v x > 0 
if E G {[ 1 / 7 $], A x , A y , A z } 


(40) 


As for the positive/negative y and z directions, the procedure is the same, replacing 
v x v y and d 7 » z , respectively. 

In this way, linear extrapolation outer boundary conditions are applied to the vector 
potential variables {[•^/ 7 $], A*} and zero-derivative, outflow outer boundary conditions 
are applied to the hydrodynamic variables {P, Po,v 1 }. These conditions are applied 
to the innermost gridpoints on the three-gridpoint-thick outer boundary surface, first 
in the positive x, then y, then z directions, followed by the negative x, then y, then 
z directions. Next, they are applied to the second innermost gridpoints on the outer 
boundary surface in all directions, and finally to the outermost point. 

Currently, IllinoisGRMHD supports the use of the xy— plane as a symmetry plane, 
in which case the negative z—direction outer boundary condition is not applied, letting 
the Cactus/Carpet parallel AMR infrastructure impose the reflection symmetry. 


3.6. Linkage of IllinoisGRMHD to the Rest of the Einstein Toolkit 

In order to evolve the GRMHD equations in a dynamical spacetime context, 
IllinoisGRMHD must be coupled to a separate module that evolves the spacetime metric, 
typically using components of the stress-energy tensor produced by IllinoisGRMHD as 
source terms. The ET is based within the Cactus infrastructure, thus modules are 
called “thorns”, of which IllinoisGRMHD is one. ET is structured so that thorns 
evaluating the spacetime metric evolution equations (i.e., the left-hand side of Einstein’s 
equations) must couple to a common interface thorn, called ADMBase. Similarly, thorns 
that evaluate evolution equations governing the right-hand side of Einstein’s equations 
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couple to an interface thorn called TmunuBase. TmunuBase and ADMBase are designed 
to interface seamlessly, so that GRMHD evolution thorns coupled to TmunuBase will 
automatically work with any spacetime evolution thorn properly coupled to ADMBase. 

Thus, since IllinoisGRMHD is fully coupled to TmunuBase, it is immediately compatible 
with all spacetime evolution formulations within ET, including BSSN, and conformal 
and covariant Z4 |lj (both provided by the McLachlan [15, US] Thorn). For dynamical 
spacetime evolutions within this paper, IllinoisGRMHD is coupled to the McLachlan 
BSSN thorn. 

4. Code Validation Tests 

This section compares the results of IllinoisGRMHD to those of two other codes written 
using the ET infrastructure: GRHydro, which is the only other open-source GRMHD 
code within ET, and the original, closed-source GRMHD code of the Illinois group 
(OrigGRMHD), on which IllinoisGRMHD is based. OrigGRMHD has been subjected 
to a large battery of stringent test-bed problems, including but not limited to standard 
ID relativistic MHD shock tests, 2D cylindrical blast explosion tests, magnetized Bondi 
accretion and stellar collapse tests as well as (self-)convergence tests [30] . Though 
it may be argued that IllinoisGRMHD has not been as robustly tested as GRHydro 
or OrigGRMHD, we demonstrate here that IllinoisGRMHD and OrigGRMHD results 
agree to roundoff error, indicating that both are algorithmically identical, and GRHydro 
and IllinoisGRMHD results agree within truncation error, indicating that both can be 
expected to converge to the same result. 

The first two parts of this section (Secs. 14.11 and 14.21) demonstrate that 
IllinoisGRMHD and OrigGRMHD generate output identical to roundoff error, in two 
complementary, highly-challenging tests. In the first test (Sec. l4.ll) . a weakly-magnetized 
TOV star is evolved over many dynamical timescales and grid light-crossing times. The 
second test (Sec. 14.2ft exposes both codes to a type of “fuzzing”, in which random 
initial data are evolved. Unlike the first test, initial data in the second test contain 
strong shocks, highly-magnetized and highly-relativistic, stochastic flows, as well as 
nontrivial, discontinuous spacetime-metric and extrinsic curvature components. Despite 
the harshness of the second test, both codes are shown in Sec. 14.21 to produce roundoff- 
error identical results over many grid light-crossing times. Results from these tests are 
highly significant, as they demonstrate that IllinoisGRMHD yields identical results to 
the “battle-hardened”, trusted code on which it is based. 

In our final validation test, both IllinoisGRMHD and GRHydro evolve 
unmagnetized, stable TOV stars in a dynamical spacetime backdrop, and are shown 
to converge to the same result at the expected order, though IllinoisGRMHD exhibits 
slightly slower central density drift and lower Hamiltonian constraint violations at a 
given resolution. 
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4-1. IllinoisGRMHD and OrigGRMHD: Roundoff-Error Agreement in Evolving 
Magnetized TOV Star 


IllinoisGRMHD, OrigGRMHD, and GRHydro all implement double-precision, 64-bit 
floating point arithmetic, which represents numbers to between 15-17 significant digits. 
Given the sophisticated and iterative nature of these GRMHD codes, initial machine- 
precision differences can grow enormously over time. As an example, we multiply the 
initial rest-mass density of a weakly-magnetized TOV star by 1 + 10“ 15 , yielding a 15th 
significant digit perturbation. We then perform the evolution, measuring the number of 
significant digits of agreement between this perturbed run and an unperturbed evolution, 
through 15 dynamical timescales and on AMR grids with multiple levels of refinement. 
The dashed blue line of Fig. [j] plots the result from this test. Notice that the number 
of significant digits of agreement quickly drops from 14 digits, plateauing to between 6 
and 8 digits of agreement. 

We were careful to develop IllinoisGRMHD so that its results agree with 
OrigGRMHD to roundoff error, and Fig. |T| confirms that for this weakly-magnetized 
TOV star test, the number of significant digits of agreement between IllinoisGRMHD 
and OrigGRMHD (solid red line) follow the same curve as the expected roundoff error 
intrinsic to OrigGRMHD (dashed blue line). For full details of the physical scenario 


modeled here, as well as the grid parameters, see Appendix B 


4-2. IllinoisGRMHD and OrigGRMHD: Roundoff-Error Agreement in Evolving 
Random Initial Data 

Although we have demonstrated that when evolving weakly-magnetized TOV stars on 
AMR grids, IllinoisGRMHD and OrigGRMHD produce results that agree to roundoff 
error, one might argue that even though there is a sharp discontinuity at the stellar 
surface, this code test is insufficient for truly demonstrating roundoff-level agreement, as 
it lacks strong shocks and highly-relativistic, highly-magnetized fluid flows. To address 
this potential criticism, we developed a random initial data module that sets up both 
weak and strong stochastic GRMHD flows atop an artificially-static, weak and strong- 
field stochastic spacetime background that is nearly conformally flat. The stochastic 
nature of these data means that both metric and GRMHD quantities suffer from both 
weak and strong discontinuities from one spatial point to the next, providing a robust 
test of the high-resolution shock-capturing algorithms within these GRMHD codes, as 
well as a confirmation of the stability of both GRMHD codes to fuzz testing. We stress 
that although the chosen metric has Lorentzian signature and the spatial three-metric 
is positive-definite, these initial data are for numerical convenience only and are not 
designed to satisfy Einstein’s equations. In fact, when these GRMHD data are evolved 
forward in time, spacetime field variables are strictly held fixed in time. 

All components of spacetime and GRMHD tensors and vectors are nonzero, 
randomly fluctuating from one spatial point to the next, with each component having 
a unique magnitude. As a result, this module has been useful in checking for typos in 
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Roundoff-Error Code Agreement: 
Dynamic Spacetime GRMHD TOY Test 



Figure 1 . Significant digits of agreement between pairs of codes, monitoring 
the central density of a magnetized neutron star in a dynamic-spacetime GRMHD 
simulation versus time, as measured in dynamical timescales fdyn = l/^/po,max- 
The solid red line shows the number of significant digits of agreement between 
IllinoisGRMHD and OrigGRMHD. The dashed blue line shows the expected roundoff 
error, measured as the number of significant digits of agreement between OrigGRMHD 
and itself with a 15th-significant-digit perturbation to the initial density of the 
magnetized neutron star. This run was performed on 8 parallel processes on a desktop 
computer. All details from this simulation are provided in|Appendix B| 


the GRMHD evolution equations, which were completely rewritten in IllinoisGRMHD. 
For example, if by mistake 7 xy were written 7 xz in any of the GRMHD equations, then 
because these components have differing magnitudes, IllinoisGRMHD and the original 
GRMHD code of the Illinois group would not agree to roundoff precision and the test 
would fail. This module was used extensively in the first stages of IllinoisGRMHD 
development to find such typos, as well as truncation-error-level algorithmic differences 
between the old and new codes. All such typos and algorithmic differences were fixed and 
modified, respectively, so that roundoff-level agreement could be demonstrated. A full 


description of the random initial data module and grid setup is provided in Appendix C 


As with the magnetized TOV roundoff-error test of Sec. 14.11 we measure the 
expected level of roundoff error by first adding a random, 15th-digit perturbation to 
all GRMHD primitive variables after they are set. Then we evolve both perturbed and 
unperturbed initial data with the trusted OrigGRMHD code. The difference grows with 
time, but plateaus to about 13 digits over time, as shown in Fig. [2j 

Next, the same unperturbed initial data are evolved on the same grids with 
IllinoisGRMHD, and the results confirm that IllinoisGRMHD and OrigGRMHD indeed 
agree to within expected roundoff error, through at least 25 light-crossing times. All 
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Roundoff-Error Code Agreement: 
Random Initial Data 



Figure 2. Significant digits of agreement between pairs of codes, monitoring the rest- 
mass density po summed along the x-axis on the finest AMR level. The solid red 
line shows the number of significant digits of agreement between IllinoisGRMHD and 
OrigGRMHD. The dashed blue line shows the expected roundoff error, measured as 
the number of significant digits of agreement between OrigGRMHD and itself with a 
random 15th-significant-digit perturbation to all primitive GRMHD variables. This 
run was performed on 8 parallel processes on a desktop computer, and was shown to 
agree with the single-process OrigGRMHD run to roundoff-error as well. All details 
from this simulation are provided in |Appendix B| 


of these runs were performed on 8 parallel processes. Full details regarding how these 
initial data are generated, as well as the computational setup for these simulations, are 
provided in Appendix B 


4.3. IllinoisGRMHD and GRHydro: Unmagnetized TOV Star Convergence Tests 

The open-source IllinoisGRMHD and closed-source OrigGRMHD have been shown to 
produce roundoff-error identical results even when evolving very harsh, relativistic, 
strongly-magnetized, discontinuous initial data. IllinoisGRMHD represents the second 
open-source, dynamical spacetime GRMHD module, the first being GRHydro [50]. Both 
are based in the Einstein Toolkit, which provides a particularly convenient infrastructure 
for performing GRMHD simulations in a dynamical spacetime context. GRHydro 
contains a large number of features, including a variety of reconstruction options, 
approximate Riemann solvers, and outer boundary options. OrigGRMHD contains 
many such features as well, but nearly all of these features are not robust in the 
context of black-hole-inhabited spacetimes and have thus remained unused for years. 
IllinoisGRMHD contains only the features from OrigGRMHD that have been used in 
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all recent papers by the Illinois NR group (e.g., [2SJ [-ST. SHI El ISO]). 

This section compares results between IllinoisGRMHD and the standard, 
FORTRAN version of GRHydro, using identical initial data, computational grids, 
dynamical spacetime evolution (BSSN) modules, reconstruction scheme, and Riemann 
solver. Since IllinoisGRMHD has been shown to agree with OrigGRMHD to roundoff 
precision, these tests can also be seen as a proxy comparison between GRHydro and 
OrigGRMHD. 

As detailed in Appendix A both IllinoisGRMHD and GRHydro evolve the same 
physical scenario in this test as in Sec. 14.11 but with the magnetic fields inside the 
TOV star set to zero. GRHydro options were chosen so that its algorithms would be 
identical to IllinoisGRMHD. Despite the basic algorithms being the same, both codes 
differ significantly in how they are implemented. This difference in implementation 
should result in truncation-level differences between the two codes, but instead we find 
slightly different convergence properties between the codes. 

A quantity Q(Ax) that converges to zero at nth order with increasing resolution 
(i.e., decreasing grid spacing Ax) satisfies 


Q( Asi) _ /AxA" 

Q(Ax 2 ) \Ax 2 J 

Thus the convergence order to zero, n, is written as follows: 


f Q(Ax i) \ 


" = l0g U(A^J/ log 


Axi 

Axo 


(42) 


Figure [3] demonstrates that for this stable, equilibrium TOV star, truncation 
errors lead to nonzero drifts in star’s central density and the L2-Norm of the 
Hamiltonian constraint, which each converge to zero at roughly second-order [n(t) ~ 2, 
where n is as defined in Eq. (1421) ]. Notice that L2-Norm Hamiltonian constraint 
convergence order fluctuates significantly in GRHydro evolutions, as compared to 
IllinoisGRMHD. Additionally, at the highest resolution chosen (resolving the NS 
diameter to approximately 80 gridpoints), GRHydro possesses roughly 8% higher 
Hamiltonian constraint violation than IllinoisGRMHD, as shown in Fig. [4] We conclude 
that IllinoisGRMHD appears to suffer from less Hamiltonian constraint violation than 
GRHydro at a given resolution, and exhibits more consistent L2-Norm constraint 
violation convergence to zero as resolution is increased. 

However, at all resolutions, the absolute value of central density drift through 55 
dynamical timescales is far higher with IllinoisGRMHD than GRHydro. We analyze the 
drift at high resolution in the lower panel of Fig. [4J finding that the differences appear 
directly after initial settling of the TOV star. These differences are not surprising, 
given the unique algorithmic choices in each code (e.g., GRHydro adopts the internal 
energy e as a primitive variable, where IllinoisGRMHD adopts pressure P instead, just 
to name one). In this plot, we fit data in the range 5 < t/Oyn < 55 to a least-squares 
linear trendline, finding that the rate of central density drift (i.e., the slope of the linear 
trendline) after 5 dynamical timescales to be within about one standard deviation for 
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IllinoisGRMHD: 

Ap c (t) Convergence to Zero 




IllinoisGRMHD: 



Figure 3. IllinoisGRMHD and GRHydro convergence tests, for dynamic-spacetime, 
unmagnetized equilibrium TOV star evolutions (with physical and numerical setup 
as described in | Appendix A[ ). Upper panels: Convergence to zero of TOV star 
central density drift A p c (t) = p c (t)/p c ( 0) — 1, comparing GRHydro (left plot) with 
IllinoisGRMHD (right plot). Top plots show A p c (t) at three separate resolutions, 
with the low (dotted magenta) and medium (dashed blue) resolution (LR and MR, 
respectively) simulation results rescaled to high resolution (HR, solid red), assuming 
that A p c {t) converges to zero at second order. Lower plots show implied convergence 
order to zero [see Eq. (l42l) ] of A p c (t) for pairs of runs, for HR and MR (thin dashed 
black), and MR and LR (thick solid red), where convergence order to zero is defined as 
in Eq. (l42l) . Lower panels: Convergence to zero Eq. (l42l) of L2 Norm of Hamiltonian 
constraint violation, ||iJ(i)||. Top plots show ||iL(t)|| at three resolutions, rescaled so 
that LR (dotted magenta) and MR (dashed blue) results should overlap HR (solid red) 
results if second-order convergence to zero is achieved. The bottom plots show implied 
observed convergence order to zero of pairs of runs: HR and MR (dashed black), and 
MR and LR (dashed blue). 


the two codes. In addition, we verified that although the simulation is run for about 
2.3 light-crossing times, doubling the outer boundary has no qualitative effect on the 
results. 

In a forthcoming paper, we will demonstrate that differences between these two 
open-source GRMHD codes spawn from how the GRMHD evolution algorithms are 
implemented, independent of the chosen reconstruction scheme. 
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Figure 4. Truncation-error analysis, comparing results from IllinoisGRMHD and 
GRHydro at high resolution (HR). The top plot shows L2-Norm Hamiltonian 
constraint violation, for IllinoisGRMHD (red solid) and GRHydro (blue dashed). 
GRHydro exhibits about 8% higher constraint violation, so its data were multiplied 
by 0.92 to achieve a good overlap with IllinoisGRMHD data. The bottom plot shows 
central density drift A p c (t) = p c (t)/p c { 0) — 1, at high resolution (HR) as well, for 
IllinoisGRMHD (thin dashed blue) and GRHydro (thin solid red). The thick blue 
and red lines are linear least-squares fits to IllinoisGRMHD and GRHydro data, 
from 5tdyn to 55tdyn- The slope on the GRHydro line is (—1.19 ± 0.02) x 10 -5 , 
and (—1.13 ± 0.04) x 10 -5 for IllinoisGRMHD, where the errors given are standard 
deviations. 


5. Performance Benchmarks 

We have demonstrated that although IllinoisGRMHD represents a complete rewrite of 
OrigGRMHD, the two codes agree to roundoff precision. IllinoisGRMHD is also designed 
to be more user-friendly, more extensible, and better documented than OrigGRMHD 
as well. Furthermore, as we demonstrate in this section, IllinoisGRMHD performs 
and scales better than OrigGRMHD. This stems from the fact that coding decisions 
within IllinoisGRMHD were made specifically from the outset to optimize not only 
user-friendliness and code readability, but also performance. 

Making IllinoisGRMHD perform as well as GRHydro, on the other hand, 
appears to be an unlikely goal, as the AMR-capable GRMHD algorithm adopted by 
IllinoisGRMHD/OrigGRMHD is far more computationally intensive. All variables in 
GRHydro’s GRMHD scheme for AMR grids (hyperbolic divergence cleaning I l9l [50] ) 
are unstaggered, overlapping gridpoints. Meanwhile, IllinoisGRMHD/OrigGRMHD 
implement a staggered vector-potential formulation, requiring, e.g., about 60% more 
expensive reconstructions to compute GRMHD fluxes, as they must be computed 
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on staggered gridpoints. In addition, the evolution of the staggered EM vector 
potential gauge quantity [^^7$] is quite expensive, as it requires a large number of 
interpolations. Of course, in exchange for this more expensive algorithm is the guarantee 
that monopoles (i.e., violations of V • B = 0) cannot be generated on grid refinement 
boundaries when magnetized fluid flows cross them. GRHydro cannot guarantee this, 
but IllinoisGRMHD/OrigGRMHD does. 

Thus a priori , we would expect IllinoisGRMHD and OrigGRMHD to significantly 
under-perform GRHydro. Remarkably, Fig. 0 demonstrates that for a physical system 
and AMR grid hierarchy typically used in production runs, IllinoisGRMHD actually 
outperforms both GRHydro and OrigGRMHD by a comfortable margin. There 
also exists a new, experimental C++ version of GRHydro (henceforth, GRHydro- 
experimental), which was written in part to improve performance. Indeed, performance 
is greatly enhanced by GRHydro-experimental, but it at best matches IllinoisGRMHD 
performance at small core coimts, and scales worse than IllinoisGRMHD with increasing 
core count. The physical system and basic AMR grid hierarchy is as described in 
Appendix A[ i.e., it consists of an unmagnetized, equilibrium TOY star for which the 


magnetic-field is also evolved but initialized to zero. 

As measured by the number of gridpoints computed per second per core on the 
TACC Stampede supercomputer, at all problem scales typically used for parallel AMR 
runs (ranging from 32 to 2,048 cores), Fig. [5] shows that IllinoisGRMHD consistently 
outperforms the standard GRHydro by a factor of between 1.7-1.8. However, 
IllinoisGRMHD matches GRHydro-experimentars performance, to within measurement 
error at 32 cores, but manages to outperform GRHydro-experimental by about 16% at 
2,048 cores. Again, it is remarkable that IllinoisGRMHD can produce performance 
numbers in the same ballpark as GRHydro, as IllinoisGRMHD implements a much 
more expensive GRMHD algorithm. 

The performance improvement over OrigGRMHD is also significant, with 
IllinoisGRMHD outperforming OrigGRMHD by a factor of 1.3 at 32 cores, increasing 
to 1.6 at 256 and 2,048 cores. In independent testing, we find that about 10 — 20% 
of the performance difference between IllinoisGRMHD and OrigGRMHD is due to the 
fact that OrigGRMHD is based on an old, unmaintained version of the Cactus/Carpet 
infrastructure (ca. October, 2010). For more details on the physical system and basic 


grid structure of this benchmark, see Appendix A 


Benchmarks presented here measure total simulation performance, and since these 
are dynamical spacetime simulations, the performance gap between IllinoisGRMHD 
and the other codes will certainly increase for fixed-background-spacetime simulations. 
Thus by adopting IllinoisGRMHD, research groups currently using the standard version 
of GRHydro or OrigGRMHD stand to boost their computational resources by a factor 
of between 1.6-1.8. Independent tests indicate that the performance gap increases to a 
factor of & 2 in /Ae+spacetime-background GRMHD simulations. 

Making AMR-based codes like IllinoisGRMHD, OrigGRMHD, and GRHydro scale 
well is an intrinsically difficult task. AMR greatly reduces the memory and processor 
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Relative Performance Benchmark 


I 0.8 

fa 

g 0-6 

a 

| 0.4 

«S 

0.2 
0 

10 100 1000 

Number of Cores 

Figure 5. Relative performance of IllinoisGRMHD (blue solid), OrigGRMHD (red 
dashed), standard GRHydro (green dotted), and the experimental C++ version 
of GRHydro (magenta dot-dashed) at multiple problem scales on the Stampede 
supercluster, performing an unmagnetized neutron star simulation, but with magnetic 
field evolution enabled. Performance data are normalized to IllinoisGRMHD two-node 
(32-core) performance (as measured by the number of gridzones computed per second 
per processor core). As the number of cores was increased, the number of gridpoints 
per core was kept fixed at approximately 72 3 for all four refined levels and 68 3 for the 
lowest-resolution level on these AMR grids, effectively making this a weak-scaling test. 



GRHydro + 
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overhead in our simulations, focusing resolution only where it is needed, and generating 
many small, refined numerical grids in the process. When these small refined grids 
are parallelized, however, they are generally split into even smaller grids, resulting in 
a large grid surface area to volume ratio. As the information on the surfaces must be 
communicated across nodes, this makes the performance of AMR-based codes strongly 
network-limited. 

OpenMP [53] can be used to combat this by splitting computational loops over 
multiple processor cores, enabling us to use fewer parallel (MPI) processes per CPU and 
thus larger grids on a given (MPI) process. This reduces the network load significantly 
and thus increases overall performance. 

As shown in Fig. [6j in production-scale benchmarks, we find that IllinoisGRMHD 
performs slightly more than 40% faster as an OpenMP/MPI hybrid code than as a pure 
MPI code (i.e., when running 16 MPI processes per node, OpenMP disabled), for a 
typical GRMHD production run on the Stampede supercluster. Although all three codes 
possess some degree of OpenMP support, all of IllinoisGRMHD’s loops have has been 
written with full OpenMP [53] support, making IllinoisGRMHD a pure OpenMP/MPI 
hybrid code, just like OrigGRMHD. We finish this section by noting that all benchmark 
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Benefits of OpenMP/MPI Hybridization, 
1,024-Core Test on Stampede 



Figure 6. IllinoisGRMHD code performance as the number of MPI processes per 
node on Stampede is varied, with the total core-count fixed at 1,024 (i.e., 64 Stampede 
nodes). When running with 1,2,4, 8 , and 16 MPI processes per node, the number of 
threads per MPI process (OMP_NUM_THREADS) was set to 16,8,4,2, and 1, respectively. 
Efficiency is normalized to the 4 MPI processes per node case. Simulation is of a 
neutron star with full GRMHD and spacetime evolution enabled. 


results in Fig. [5]were performed using 4 MPI processes per node and 4 OpenMP threads 
per MPI process, which we found maximized performance in all codes. 

6. Conclusions and Future Work 

The held of numerical relativity has matured considerably in the years since the first 
dynamical spacetime GRMHD codes were developed, and multiple groups now possess 
such codes. Given that the future of our held depends on the ability to advance 
and extend these codes to model new physics, while still maintaining and improving 
the GRMHD modules, it stands to reason that the community could benefit if we 
consolidated our efforts and adopted the same dynamic-spacetime GRMHD code. 

With its proven robustness and reliability in modeling some of the most extreme 
phenomena in the Universe, it seems the OrigGRMHD code could be a good candidate 
for such community adoption if it were open-sourced. But despite its strong scientihc 
track record, OrigGRMHD was not written with community adoption in mind, instead 
being a code written “by experts and for experts” of OrigGRMHD, with a premium put 
on immediate applications. As such, the code lacked a number of features common to 
large, open-source, community-based codes in computational astrophysics, including 
sufficient documentation and code comments, hne-grained modularity, a consistent 
coding style, and regular, enforced code maintenance (e.g., removal of unused and 
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unmaintained features). 

As an open-source, from-the-ground-up rewrite of OrigGRMHD, fllinoisGRMHD 
aims to fix all the former code’s idiosyncracies, thus facilitating widespread community 
adoption. With such adoption in mind, fllinoisGRMHD’s development has been 
guided by the four core design principles of user-friendliness, modularity/extensibility, 
robustness, and performance/scalability. Regarding user-friendliness, the code is well- 
documented, properly commented, and requires only basic programming skills to 
understand and run. fllinoisGRMHD is also far more modular and extensible than 
OrigGRMHD, with low-level CFD routines split off from the main code into a library 
of extensible functions. 

As for robustness, fllinoisGRMHD was designed to act as a drop-in replacement for 
OrigGRMHD, and we have demonstrated that fllinoisGRMHD indeed reproduces results 
from the original code to roundoff-level precision, not only when evolving magnetized 
neutron stars, but also discontinuous, random initial data. 

In addition, fllinoisGRMHD largely produces consistent results with the only 
other open-source, dynamical spacetime GRMHD code, GRHydro, in that both codes 
exhibit approximate second-order convergence. Although both codes were run with the 
same basic evolution algorithms, results differ due to the specific details of how these 
algorithms were implemented. We will explore this further in a forthcoming work, but 
just to name a couple of differences, GRHydro reconstructs the specific internal energy e 
and the Valencia-formulation 3-velocity, while fllinoisGRMHD reconstructs pressure and 
3-velocity defined as v l = u l /u°. When evolving equilibrium unmagnetized TOV stars, 
GRHydro produces about 8% higher Hamiltonian constraint violations (as measured 
by the L2 Norm over the entire grid), but significantly less absolute central density 
drift than fllinoisGRMHD. We find that the rate of central density drift is identical 
between the two codes after the star undergoes an initial settling over a few dynamical 
timescales. We conclude it is this initial settling that causes the large discrepancy in 
absolute central density drift. 

Though user-friendliness, modularity/extensibility, and robustness were the 
primary considerations in HlinoisGRMHD’s development, it would be hard to convince 
key developers of other codes to adopt fllinoisGRMHD unless we could demonstrate 
at least comparable performance and scalability to alternative dynamical spacetime 
GRMHD codes. To this end, we have shown that fllinoisGRMHD is in fact about 
1.7-1.8 times faster than the standard version of GRHydro for production-size, AMR- 
enabled, GRMHD runs on the Stampede supercluster, scaling to typical high-resolution 
core-counts at better than 95% efficiency. 

This is a rather remarkable result, as fllinoisGRMHD implements a far more 
computationally expensive GRMHD algorithm than GRHydro to ensure the no¬ 
monopole constraint V • B = 0 is satisfied. This added expense forbids the generation 
of monopoles at AMR grid boundaries in the case of multi-scale GRMHD flows, 
which GRHydro, and even its new, experimental C++ version, which we refer to as 
GRHydro-experimental, cannot guarantee. GRHydro-experimental, which was actively 
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being written during the preparation of this paper, is a complete rewrite of the 
standard, FORTRAN-based GRHydro, with a goal in part of improving performance. 

Indeed, GRHydro-experimental does improve performance, but only at best matches 
IllinoisGRMHD’s performance at small core counts, with IllinoisGRMHD performing 
about 16% better at 1,024 to 2,048 cores. Though these results may change when 
simulating other systems of interest, we consider them representative. 

IllinoisGRMHD also outperforms OrigGRMHD by a factor of 1.3-1.6. Thus by 
adopting IllinoisGRMHD, research groups using GRHydro or OrigGRMHD stand to 
increase their computational resources available for dynamical spacetime, GRMHD runs. 

While in terms of performance, IllinoisGRMHD seems to be only slightly better than 
the experimental, C++ version of GRHydro, perhaps the two greatest advantages of 
IllinoisGRMHD over GRHydro is that (1) IllinoisGRMHD does not allow the generation 
of magnetic monopoles when modeling multi-scale GRMHD flows on AMR grids and 
(2) IllinoisGRMHD is capable of stably modeling GRMHD flows into BH horizons over 
very long timescales, without the need for special algorithms that excise GRMHD data 
with the BH. GRHydro requires excision to model such flows and its GRMHD features 
have been mostly used for core collapse (to a neutron star) simulations, in which no 
BH is present. We conclude that making GRHydro’s GRMHD schemes as robust may 
require careful specification of boundary conditions on the excision surface coupled to 
an interpolation scheme across AMR level boundaries that respects the no-monopoles 
constraint. 

As mentioned previously, a forthcoming paper will analyze how differences in 
algorithmic implementations between GRHydro and IllinoisGRMHD can lead to 
significantly different results when evolving neutron stars. One possibility is that we may 
find an implementation that results in a superior code to either original code. If such 
a code is found, it may prove quite useful to reliably evolving binary neutron stars and 
black hole-neutron stars over many orbits with a minimum of computational expense. 

Although IllinoisGRMHD is ready for production runs now, we encourage other 
developers to join our effort in improving IllinoisGRMHD beyond its current state, as a 
great deal of important work remains to be done. We would like to port features from 
GRHydro into IllinoisGRMHD’s library of functions, using IllinoisGRMHD’s coding 
style, including reconstruction schemes, conservative-to-primitives solvers, and more 
advanced approximate Riemann solvers, just to name a few. IllinoisGRMHD was 
originally written in a standalone sandbox for maximum portability, and was only 
recently ported into the Einstein Toolkit. It therefore makes minimal use of certain 
aspects of the ET infrastructure that could greatly extend its usefulness. For example, 
the current version supports only single gamma-law EOSs, and full 3D simulations with 
either no symmetries enabled, or simply bitant symmetry across the xy- plane. The ET 
infrastructure provides support for arbitrary EOSs, symmetry conditions, etc., and we 
intend to work with the large ET community toward making this code the standard 
choice in the ET, and one the community can be proud of. 
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Appendix A. Physical and Computational Setup for Unmagnetized TOV 
Star Convergence Tests 

IllinoisGRMHD and GRHydro evolve unmagnetized, polytropic TOV star initial data, 
consisting of a TOV star with central density of 0.129285309 constructed with a 
polytropic equation of state P = Kp r , where K = 1 and T = 2. This generates a 
model for a cold, degenerate neutron star (NS) with compaction M ns /Rns ~ 0.1467. 

The initial data are generated using the built-in TOVSolver thorn within the ET. 

The convergence tests are dynamical spacetime tests, in which IllinoisGRMHD and 
GRHydro are each coupled to the McLachlan BSSN thorn. To match lllinoisGRMHD’s 
evolution algorithms, the “HLLE” approximate Riemann solver is chosen for GRHydro 
evolutions, along with PPM reconstruction. Though the codes differ in their usage of 
these algorithms (e.g., GRHydro reconstructs the internal energy e instead of pressure 
P, etc.), the parameters are chosen so that the both codes share precisely the same 
algorithms. Note that the outer boundary conditions for the spacetime variables are set 
to be identical between the two codes, but the hydrodynamic boundary conditions differ. 
However, this is inconsequential, as the density near the outer boundary remains within 
~ 1% of the original atmosphere value in both codes throughout the entire evolution. 

These tests are performed on cubic AMR grids, so that the coarsest grid cube 
possesses a half-side-length of 40.Rns, centered on the NS. The AMR hierarchy—nested 
within this coarse grid and also centered on the NS—consists of four, progressively 
higher-resolution cubes with half-side-lengths of 15RnS) 7.5Rns, 3.75Rns> and 1.875-Rns- 
At each finer level, the grid spacing is halved, so that the cube with half-side-length 

1.875-Rns has a grid spacing of Ax = {0.03906, 0.03125, 0.025}Rns for low, medium, 

and high-resolution runs, respectively. This corresponds to resolving the NS across 
its diameter to approximately 51, 64, and 80 gridpoints, for low, medium, and high- 
resolution runs, respectively. In both codes, low-density atmosphere density floor is set 
to correspond to 10”' times the initial central density of the NS. 
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Appendix B. Physical and Computational Setup for Magnetized TOV Star 
Code Validation Tests 


The magnetized TOV star used in these tests is precisely the same stable neutron star 


as described in Appendix A except it is seeded with a weak magnetic field at t — 0. 
This initial magnetic field is purely poloidal, with vector potential components defined 


as 


A x = - yA b max(P - P cut , 0), (B.l) 

A y = xA b max(P - P cut , 0), (B.2) 

A z , [«/y$] = 0, (B.3) 

where P cut is set to 4% of the maximum pressure and A b is set so that the maximum 
initial magnetic-to-gas pressure ratio /R 1 is 0.83 x 10 -3 . Like the convergence tests, these 
tests are dynamical spacetime tests, but unlike the convergence tests, IllinoisGRMHD 
and OrigGRMHD codes are each coupled to the (closed-source) BSSN dynamical 
spacetime module of the Illinois NR group. Further, the initial TOV star density, 
pressure, and spacetime metric profiles are provided by the code of G9, which generates 
data for this nonrotating star that agrees to machine precision with the TOVSolver code 
used in the unmagnetized tests. 

These tests are also performed on AMR grids, with the coarsest grid cube extended 
so its half-side-length is 10Pns> centered on the NS. Two finer AMR grid levels are 
nested within this coarse grid: the next finer grid having half-side-length 5Pns and the 
finest having 2.5 Pns- Again, at each finer level, the grid spacing is halved, so that 
the cube with half-side-length 2.5Pns has a grid spacing of Ax ~ {0.078}Pns- This 
corresponds to resolving the NS to 26 gridpoints across its diameter. The resolution is 
intentionally set to be very low, to guarantee that truncation-error differences between 
the codes will be more strongly magnified than at (higher) resolutions typically chosen 
for evolving NSs. In addition, a close outer boundary is chosen to maximize its influence 
on the evolution, to check for errors in coding the outflow outer boundary conditions. 
Finally, we choose our low-density atmosphere density floor to correspond to ICR' times 
the initial central density of the NS. 


Appendix C. Random Initial Data 

The algorithm used within our random initial data module is as follows. At a given 
gridpoint, the random number generator is first seeded with a unique integer based on 
the coordinate of the gridpoint. The seed is then used to generate one double-precision 
random number £ E (—0.1, 0.1). Because this random number is based entirely on the 
coordinate of a gridpoint, the initial data will be consistent when evolving on identical 
grids, regardless of how the global grid is split among processors when generating and 
evolving initial data in a parallel run. 
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Table Cl. Prescription for setting GRMHD and spacetime metric quantities in 
random initial data module. Note that the initial data employ the polytropic EOS 
P = Pq, and we evolve with the gamma-law EOS P = (T — l)poe with T = 2. 


Variable 

Value set in random initial data module 

Po 

0.015 

P 

Po (Gamma-law EOS) 

<t> 

0.15 

a 

1-0.15 

{v x , v y , v z } 

{1,1.1,0.9} x 10- 4 £ 

{[y/V^]) A x , Ay, A z } 

{0.6,1,1.1,0.9} x 10 _1 £ 

{fd x ,fd y ,fd z } 

{0.9,1.1,1} x 10 _1 £ 

{fxxi fyyi fzz\ 

1 +{ 10 , 1 , 10 } x 10 _3 £ 

{ 7 xyi 7 xzi fyz) 

{ 1 , 10 , 1 } x 10 _4 £ 

\K K K K K K \ 

\ 1Y xxi JV xy) J *-xz7 lv yyi 1Y yz) J ^zz f 

{10,2,3,40,5,60} x 10“ 3 £ 


Using this double-precision random number £ G (—0.1, 0.1), we define 

5 = £ + - [4 + sin( 10 x/L) + sin(10p/L) + cos( 102 r/L)], (C.l) 

where L denotes the coarsest AMR grid cube half-side-length. Thus 5 G (^, 1.1) 
consists of a strong stochastic perturbation atop a smoothly oscillating function. Given 
the quantities 5 and £, all basic spacetime metric and GRMHD quantities may then be 
set following the prescription detailed in Table IC11 

Although the resulting data are not designed to satisfy Einstein’s equations, be 
conformally flat, or be consistent with the BSSN formalism, we make sure the spacetime 
metric is Lorentzian and the three-metric positive definite. We also enforce the BSSN 
constraint 7 = 1 as follows. Immediately after the quantities in Table IG1I have been 
set at a given gridpoint, we compute the resulting non-unit determinant 7 and then 
multiply the components of 7^ by 7” 1 / 3 . In this way, the unit determinant of the three- 
metric is enforced. After applying this constraint, all GRMHD and spacetime metric 
quantities not specified in Table ICTI (e.g.. GRMHD conservative variables, B\ 7 % etc.) 
are directly computed from quantities in that table. 

In tests adopting this random initial data module, we choose an AMR grid with one 
refinement level centered at the origin. The coarser grid, with grid spacing 1.0, has cube 
half-side-length of 10, and the finer grid, with grid spacing 0.5, has grid half-side-length 
of 2 . Spacetime evolution modules, as well as prolongation and restriction operations 
on spacetime variables, are disabled. 

The initial maximum rest-mass density is chosen to be ~ 10% the rest-mass density 
of the TOV star used in other tests. The atmospheric density floor is set to 10” 6 , 
which is about 10" 4 times the maximum possible initial density, yielding stochastic 
fluctuations in density spanning about four orders of magnitude. Additionally, the 
magnetic-to-gas-pressure ratio at t — 0 ranges from ~ 10 -3 (gas-pressure dominated) 
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to about 10 (magnetic-pressure dominated) initially. Given that initial data at a given 
gridpoint are independently and randomly specified, we conclude that physical quantities 
at neighboring gridpoints can differ by several orders of magnitude, making this a very 
harsh test. Despite the lack of coherent GRMHD flows, throughout the evolution of these 
random initial data, the Lorentz factor limit of W = 10 is exceeded and subsequently 
enforced about 400 times. 
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